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Abstract 


A computational procedure is presented for the solution of frictional contact 
problems for aircraft tires. A Space Shuttle nose-gear tire is modeled using a two- 
dimensional laminated anisotropic shell theory which includes the effects of varia- 
tions in material and geometric parameters, transverse-shear deformation, and geo- 
metric nonlinearities. Contact conditions are incorporated into the formulation by 
using a perturbed Lagrangian approach with the fundamental unknowns consisting of 
the stress resultants , the generalized displacements , and the Lagrange multipliers 
associated with both contact and friction conditions . The contact-friction algorithm is 
based on a modified Coulomb friction law . A modified two-field, mixed-variational 
principle is used to obtain elemental arrays . This modification consists of augmenting 
the functional of that principle by two terms: the Lagrange multiplier vector associ- 
ated with normal and tangential node contact-load intensities and a regularization 
term that is quadratic in the Lagrange multiplier vector. These capabilities and com- 
putational features are incorporated into an in-house computer code. Experimental 
measurements were taken to define the response of the Space Shuttle nose-gear tire to 
inflation-pressure loads and to inflation-pressure loads combined with normal static 
loads against a rigid flat plate . These experimental results describe the meridional 
growth of the tire cross section caused by inflation loading, the static load-deflection 
characteristics of the tire , the geometry of the tire footprint under static loading con- 
ditions, and the normal and tangential load- intensity distributions in the tire footprint 
for the various static vertical loading conditions. Numerical results were obtained for 
the Space Shuttle nose-gear tire subjected to inflation pressure loads and combined 
inflation pressure and contact loads against a rigid flat plate. The experimental mea- 
surements and the numerical results are compared. 


Introduction 

Contact-friction problems are inherently nonlinear 
and path dependent. Nonlinearity occurs partly because 
both the contact area and the contact-load intensities are 
not known beforehand and vary during the loading his- 
tory. Path dependency is a result of the nonconserva- 
tive (irreversible dissipative) character of the frictional 
forces. 

A review of static contact problems presented in ref- 
erence 1, which includes a bibliography of approxi- 
mately 700 papers, points out that contact problems are 
important to thermomechanical stress analyses, fracture 
mechanics, mechanical problems involving elastic foun- 
dations, the mechanics of joints, geomechanics, and tires. 

Contact problems occupy a position of special 
importance in aircraft tire mechanics because the contact 
zone is where the forces are generated that support, 
guide, and maneuver the airplane. Distributions of con- 
tact loads and frictional forces define the moments and 
shears that are applied to the landing gear system (ref. 2). 
Under rolling conditions, the distribution of sliding 
velocities within the tire footprint combined with the 
frictional forces developed by the tire defines the rate of 
energy dissipation associated with the loading conditions 
and provides a measure of tire wear (refs. 3 and 4). In the 
case of the Space Shuttle orbiter, this wear mechanism is 


strong enough to cause tire failures during individual 
landing operations (refs. 5 and 6). Therefore, an under- 
standing of these tire friction forces and the resulting slip 
velocities is critical to the design of aircraft tires for the 
next generation of high-performance aircraft, such as the 
National Aero-Space Plane and the High-Speed Civil 
Transport. 

Modeling contact phenomena in the tire footprint is a 
formidable task partly because of difficulty of modeling 
tire response. Distribution of tractions and the footprint 
geometry are both functions of normal, frictional, and 
inflation tire loads. Moreover, the complex mechanisms 
of dynamic friction, which allow the tire to develop the 
necessary steering and braking forces for aircraft control 
during ground operations, are not fully understood 
(ref. 7). The tire analyst thus is forced to choose among 
several friction theories. When the tire contact problem 
includes frictional effects, the solution becomes path 
dependent and a unique solution is not guaranteed. 

The aircraft tire is a composite structure of rubber 
and textile constituents that exhibit anisotropic and non- 
homogeneous material properties. Normal tire operating 
conditions create loads that can produce large deforma- 
tions. Elevated operating temperatures from the com- 
bined effects of material hysteresis and frictional heating 
can cause variations in the material characteristics of the 



tire constituents (refs. 8-10). The laminated carcass of 
the aircraft tire is thick enough to allow significant 
transverse-shear deformations. 

These facts and attendant difficulties emphasize the 
need to develop modeling strategies and analysis meth- 
ods that include efficient, powerful and economic contact 
algorithms. Intense research has recently focused on non- 
linear analyses of static and dynamic problems involving 
contact. Novel techniques that have emerged from these 
efforts include semianalytic finite-element models for 
nonlinear analysis of shells of revolution (refs. 1 1 
and 12), reduced methods (refs. 13 and 14), and operator 
splitting techniques (refs. 15-17). References 14, 17, 
and 18 summarize applications of these new tire model- 
ing techniques. 

Objectives and Scope 

NASA Langley Research Center tire modeling 
research concentrates on developing an accurate and effi- 
cient strategy for predicting aircraft tire responses to a 
variety of loading conditions. This research focuses on 
developing tire contact modeling techniques, and the 
specific objectives of this research are (1) to develop a 
contact algorithm with friction effects included to predict 
tire response to combined inflation-pressure and static 
vertical-loading conditions, (2) to demonstrate the capa- 
bilities of this algorithm through numerical studies, and 
(3) to validate these numerical results with experimental 
data. Distribution of normal and frictional forces in the 
tire contact zone (or footprint area) is of particular 
interest. 

The contact algorithm is incorporated into a mixed- 
formulation, two-field, two-dimensional finite-element 
model based on the moderate-rotation Sanders- 
Budiansky shell theory, including the effects of 
transverse-shear deformations, laminated anisotropic 
material response, and nonhomogeneous shell character- 
istics (refs. 19 and 20). A perturbed Lagrangian formula- 
tion (refs. 21 and 22) is the basis for this contact 
algorithm. The Lagrangian formulation uses the precon- 
ditioned conjugate gradient (PCG) iteration procedure 
(refs. 23-25) to determine contact area, distribution of 
normal force intensities, and allocation of friction force 
intensities. A modified version of the Coulomb friction 
law is incorporated into the contact algorithm in which 
the friction coefficient at the onset of sliding differs from 
that during sliding. This algorithm also monitors the 
energy dissipated within the sliding portion of the contact 
zone. In this investigation it will be assumed that the tire 
is loaded on a surface that is much stiffer than the tire, 
thus the surface will be treated as rigid. Hence, the static 
tire contact problem will be treated as a unilateral contact 


problem. Reference 26 summarizes the characteristics of 
this algorithm. 

Numerical studies presented for an inflated Space 
Shuttle nose-gear tire under static load on a flat surface 
demonstrate the capabilities of the analysis techniques. 
These analyses incorporate both friction and frictionless 
contact. Detailed studies are made of the effects of tire 
tread pattern on the contact-force intensities, the influ- 
ence of friction coefficient variations on the distribution 
of tire contact-force intensities, the convergence charac- 
teristics of the contact algorithm, and the history of 
energy dissipation in the static footprint. 

Experimental measurements were carried out on the 
Space Shuttle orbiter nose-gear tire to define its response 
to combined inflation -pressure and static vertical-loading 
conditions. Experimental procedures used to define the 
tire structural response to loading conditions and to mea- 
sure the footprint-force intensities, and empirical proce- 
dures used to define the geometry and construction 
details of the tire for modeling purposes are discussed. 
Finally, the analytical results are compared with the 
experimental measurements. 

Reference 27 describes numerical studies, experi- 
mental measurements, and comparisons between analyti- 
cal results and experimental measurements. This report 
describes development of this contact algorithm. 


Nomenclature 



nodal contact area with variable 
weighting function (see eq. (20)) 

^node 

nodal contact area (see eq. (19)) 

^quad 

area of finite-element quadrant 
(see eq. (18) and fig. 4) 

a \\> a 22 9 a \2 9 
a ss 9 a ee> a s& 
b ip ^22’ ^12* 

shell compliance coefficients 
(see eq. (4)) 

8\ v 822'8\2 

c 

number of nodal points in contact 
within element 

c <r 

portions of shell boundary where 
tractions and displacements are 
prescribed 

C <f d V’ f‘j 

tire stiffness coefficients 

0W= 1,2,6) 

c 44’ c 45’ c 55 

transverse-shear stiffness coeffi- 
cients of tire (see eq. (A 10)) 
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^0 

tangential unit vectors in meridional 
and circumferential directions 

[F] 

flexibility matrix for an individual 
element 


vector defined in equation (11) 


vector of nonlinear terms 
(see eq. (7)) 

{ G(X )} 

vector of nonlinear terms 
(see eqs. (8) and (9)) 

(G(X)} 

vector of nonlinear contributions to 
the global equations (see eqs. (11)) 


Q s , Qq transverse- shear stress resultants 

(see fig. 1 ) 

R v R 2 principal radii of curvature in meridi 

onal and circumferential directions 

r normal distance from tire axis to ref- 

erence surface 

[5] strain-displacement matrix for an 

individual element 

s meridional coordinate of tire 

(see fig. 1 ) 

[f] transformation matrix 
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{P} 


P 0 
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[QUR] 


current gaps 
initial gaps 

vector of initial gaps for contact 
element 

vector of stress-resultant parameters 
total thickness of tire 

nondimensional thickness of tire 
(see fig. 4) 

global linear stiffness matrix 
(see eq. (11)) 

bending and twisting stress resultants 
(see fig. 1 ) 

vector of nonlinear terms 
(see eq. (7)) 

number of displacement nodes in 
element 

shape functions used for approximat- 
ing generalized displacements and 
Lagrange multipliers 

extensional stress resultants 


T 

T 


n 


r 


T T 
1 1 \) 


U 

U s\ip u ’ ^slip,, 
^slip u > ^slip v 


u, v, w 


W slip’ V slip 
w slip’ v slip 


total number of degrees of freedom 

peripheral node degrees of freedom 
(see eq. (20)) 

unit normal to reference surface 
nodal (contact) force 
load parameter 

normalized external load parameter 

global vector of normalized external 
loads and initial gaps 

intensity of inflation pressure 

intensity of external loading in coor- 
dinate directions (see fig. 1) 

elemental matrices associated with 
contact condition and regularization 
term in functional 


W 

w 


y, z 


*3 

{Z} 

A g v 

e ,’ e n 

E relax 

e 0 > 2e ,0 


intensity of contact force acting nor- 
mal to contact surface 

resultant contact friction force 
(see eq. (2a)) 

intensity of contact friction forces 
tangent to contact surface 

strain energy density (strain energy 
per unit area) 

energy dissipated during slip at each 
iteration (see eqs. (27)) 

energy dissipated during slip 
over all iterations in a load step 
(see eqs. (28)) 

displacement components of refer- 
ence surface of tire in meridional, 
circumferential, and normal 
directions (see fig. 1 ) 

tangential slip for each iteration (see 
eqs. (25)) 

total tangential slip for each load step 
(see eqs. (26)) 

external work 

vector of nodal displacements in 
shell coordinate system 

vector of nodal displacements in 
Cartesian coordinate system 

Cartesian coordinate system 

coordinate normal to tire reference 
surface (see fig. 1) 

global response vector 
slip distances in tire footprint 

penalty parameters in tangential and 
normal directions (see eq. (6)) 

relaxation parameter (see eqs. (16)) 

extensional strains of reference sur- 
face of tire 
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Superscripts: 

(O 

ij 


r 
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transverse shear strains of tire 

circumferential (hoop) coordinate of 
tire (see fig. 1) 

bending strains of tire 

Lagrange multiplier, representing 
intensity of contact load normal to 
contact surface 

Lagrange multipliers, representing 
intensity of contact friction loads 
tangent to contact surface 

Lagrange multipliers, representing 
sliding friction load intensities 

vector of nodal values of Lagrange 
multipliers 

static coefficient of friction 
dynamic coefficient of friction 
dimensionless coordinates along 
meridian (see fig. 4) 

functionals 

rotation about normal to tire refer- 
ence surface (see eq. (A9)) 

rotational components of reference 
surface of tire (see fig. 1) 

element domain 
contact surface 

angle defined by ratio of two friction 
forces (see eq. (2d)) 

= 3/95 
s 3/90 
first variation 


individual elements 

indices of shape functions for 
approximating Lagrange multipliers 

index of shape function for approxi- 
mating generalized displacements 
(i = 1 ,m) 

number of iteration cycles 
matrix transposition 


unknowns consist of five generalized displacements and 
eight stress resultants. Figure 1(a) gives sign convention 
for the generalized displacements and stress resultants 
and figure 1(b) shows a free body diagram of applied 
loads, torques, and contact forces. Fundamental equa- 
tions of the shell theory used herein are given in refer- 
ences 19 and 20 and are summarized in appendix A. 

Normal Contact Force Formulation 

Figure 2 shows the geometry of contact of a shell 
pressed against a flat surface. Figure 2(a) shows sche- 
matically normal gaps between the tire carcass and the 
flat surface. In the figure, £l c refers to the contact 
region; is the initial normal gap between the tire 
shell and the plate; g w is the current normal gap; and T ' 
is the normal traction on £l c . Both g w and g w are 
defined to be in the direction of the normal to and are 
measured relative to the inflated profile of the tire. Tire 
constraints normal to the contact surface can be 
expressed in terms of the following inequalities and 
equation that must be satisfied at each point on the con- 
tact surface : 


o 

Al 

* 

l&C 

(la) 

T n <0 

(lb) 

T nS w = 0 

(lc) 


The first inequality (eq. (la)) represents the kine- 
matic condition of no penetration of the contact surface 
(g w - 0 for the points in contact). The second inequal- 
ity (eq. (lb)) is the static condition of compressive (or 
zero) normal tractions. The third equation (eq. (lc)) 
states that there is zero work done by the normal contact 
stresses (i. e., the normal contact stresses exist only at the 
points where the tire is in contact with the rigid plate). 
The following inequalities are henceforth referred to as 
the inactive contact conditions : 


l w >0 

(Id) 

T n > 0 

(le) 


Equations (1) must be satisfied for both frictionless and 
frictional contact. 

Tangential Contact Force Formulation 


Mathematical Formulation 

The analytical formulation for contact of aircraft 
tires is based on a form of moderate-rotation Sanders- 
Budiansky shell theory and includes the effects of large 
displacements and transverse-shear deformation. A 
mixed formulation is used in which the fundamental 


Equations (1) define the normal tire contact con- 
straint conditions and are augmented to include friction 
constraints associated with the Coulomb friction law 
when friction forces are considered. The Coulomb fric- 
tion law is modified to include a static friction coefficient 
associated with the nonskidding or “stick” friction con- 
straint, and a dynamic friction coefficient associated with 
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the slip condition when the static friction constraint is 
violated. 

The “stick” condition defining the static friction con- 
straint is 


M + K)/T n <^ mic 

(2a) 

T u = cos (V)7> dynamic 

(2b) 

T v = sin(vj/)r n n dynamic 

(2c) 


represent the slip condition when equation (2a) is vio- 
lated, where T u and T v are the tangential tractions on the 
friction surface representing the friction forces. The 
angle \\f defines the ratio between the two friction force 
components for the slip condition and is expressed as: 


“ 

f y V 

arctan 

V 


f~ 

_ 



where T u and T v in equation (2d) are the computed tan- 
gential traction components that exist whenever the static 
friction constraint in equation (2a) is first violated. The 
following inequalities also hold on the friction surface: 


^dynamic — ^static 

(2c) 

T u**u* 0 

(2f) 

T v Ag v < 0 

(2g) 


Inequality (2e) simply states that the dynamic 
friction coefficient cannot be greater than the static 
friction coefficient. Inequalities (2f) and (2g) state that 
the energy dissipated by sliding is never positive since 
the friction forces always oppose slipping in the tire 
footprint. 

Figure 2(b) shows schematically the relationship 
among the various tangential gap definitions. Initial gaps 
g UQ and g V() represent the tangential displacement of the 
inflated tire from the uninflated configuration. Current 
gaps g u and g v represent the displacement of the current 
contact solution from the inflated configuration. Delta 
gaps A g u and A g v represent the tangential slip distances 
from the previous tire contact “stick” location in the tire 
footprint. For nodes not currently in contact with the flat 
plate, the initial gaps and current gaps represent the 
orthogonal projections of the gap expressions on the con- 
tact surface. Delta gaps are defined only after contact has 
been established. 


Governing Finite-Element Equations 

Discrete equations governing the response of the tire are obtained by applying a modified form of the two-field, 
Hellinger-Reissner mixed- variational principle. This principle can be expressed in the following form: 


5r W^’ *e> N& u* M sQ , Q s , fie. U, V, W, ♦e) = 8n - 8 w 


(3) 


where 
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and 


w = J ( P s u + PqV + pw) dCl + J N N s u + N sQ v + M s <b s + A/ i0 < |> e + \q s + N s [^- - d s w^ + jV J0 ^ - l - d 

c <A — 

^N s qu + N q v + M s6 $ s + Mq^q + Qq + - ^ M ') + ^eQr - r 9 0 W ') 

J^j^.s( - “ + M ) + N s q(- v + v) + M s (- ^ + <t> 5 ) + M s q(- <j > 0 + <t> 6 ) + Q s + N s ^jj- -d 


,vv 


w >n 
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dc 


(5) 


In equations (4) and (5), a, and g are shell compliance coefficients which are the inverse of the shell stiffness 
coefficients given in appendix A; <|> is the rotation about the normal to the shell and is also given in appendix A; p s , p 0 , 
and p are the intensities of the external distributed loads in the meridional, circumferential, and radial directions, respec- 
tively; £2 is the shell domain; and c a and c u are the portions of the boundary over which tractions and displacements 
are prescribed. Quantities with a tilde (~) denote prescribed boundary stress resultants and generalized displacements; 
the underlined terms in equations (4) and (5) represent nonlinear contributions; and n s and n 0 are unit normals to the 
boundary. 

Modification to the variational principle consists of augmenting the functional of that principle by two terms: the 
Lagrange multiplier associated with the nodal contact pressures and a regularization term which is quadratic in the 
Lagrange multipliers. References 21, 22, and 28 give a detailed discussion of the perturbed and the augmented 
Lagrangian formulations. 

The modified functional has the following form: 


n = n„* + Jp.£-^a) 2 ]rf£> (6) 

n c 

where is the functional of the Hellinger-Reissner variational principle; X is the Lagrange multiplier; and e is the 
penalty parameter associated with the regularization term. Note that the addition of the regularization term amounts to 
approximating the rigid plate by continuously distributed springs with stiffness of e, for sufficiently large e. As 1/e 
approaches zero, the continuous springs become the rigid plate. 

Shape functions used in approximating the generalized displacements and the Lagrange multipliers are selected to 
be the same and differ from those used in approximating the stress resultants. Moreover, because of the nature of the 
functional II in equation (6), the continuity of neither the stress resultants nor the Lagrange multiplier is imposed at the 
interelement boundaries. 


Finite-element equations for each individual element can be cast in the following compact form: 
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( 7 ) 


6 



where {//}, { X }, and {X} are the vectors of the stress-resultant parameters, nodal values of the generalized displace- 
ments, and nodal values of the Lagrange multipliers, respectively, [F] is the matrix of linear flexibility coefficients, 
[ S ] is the strain-displacement matrix, [ Q ] and [/?] are the matrices associated with the contact condition and the regu- 
larization term in the functional, respectively, {G(X)} and {M(H, X)} are vectors of nonlinear terms, and {g 0 } is the 
vector of initial gaps in the contact region Q f . A dot refers to a zero submatrix or subvector, superscript (e) refers to indi- 
vidual elements, {P} is the normalized external load vector, and p is a load parameter. As the load is incremented, only 
the value of the load parameter p changes and the normalized load vector {P} is constant. Appendix B gives the formu- 
las for the elemental arrays [P], [S], (G(X)}, {M{H, X)}, [Q], [R], {P}, and {g Q }. 

Note that the sizes of the coefficient matrices [Q], [P], and {g 0 } var y the number of active contact condi- 
tions. The difficulty associated with an equation system whose size varies during the solution process is alleviated by 
allowing the Lagrange multipliers to be discontinuous at interelement boundaries and then eliminating them on the ele- 
ment level. If the stress-resultant parameters and Lagrange multiplier parameters are eliminated from equation (7) then 
the following equations in terms of nodal displacements {X} are obtained: 


m'tFr'tsj-eieim W 


(e) 


W 


(e) + {G(X)} (< ’ ) + e[e][^r l {g 0 } (< ’ ) -p{/ , } (e) =0 


( 8 ) 


where 


{G(X)} W = [5] , [F]" 1 {G(X)} (e) + {A/(//,X)} W (9) 

and the vector { H} in {M(H, X)} is replaced by its expression in terms of {X}. Equation (8) is the tangent operator for 
the Newton-Raphson iterative solution procedure used in this investigation and is derived in detail in appendix C. 

To simplify the treatment of the contact conditions, the displacement components are transformed from shell coordi- 
nates ( s , 0, jc 3 ) to the global Cartesian coordinates ( x , y , z) before assembly. The relations between the displacement 
vector in the shell coordinates {X}^ and the corresponding vector in Cartesian coordinates {X}^ can be written in 
the following compact form: 

{X} (e) = [T]{X} ie) (10) 

where [P] is the transformation matrix. The different arrays in the finite-element equations are transformed accordingly. 
Appendix D gives the explicit form of the transformation relations. 


Solution of Nonlinear Algebraic Equations 

Discrete equations governing the response of the tire 
are obtained by assembling the elemental contributions 
in equation (6) or (7) and can be written in the following 
form: 

{/(Z, p ) } = [k}{Z} + {G(Z)}-p{~P} = 0 (11) 

where [X] is the global linear stiffness matrix of the tire; 
(G(Z)} is the vector of nonlinear contributions; {P} is 
the global vector of normalized external loads and initial 
gaps; and { Z } is the global response vector of the tire 
obtained by assembling the contributions from the sub- 
vectors { H }, { X } , and { X } . 

The nonlinear algebraic equation (eq. (8)) is solved 
and the contact region and the contact-load intensities are 
determined by using an incremental-iterative technique 
(i.e., a predictor-corrector computation method) in which 
the response vector {Z}, corresponding to a particular 
value of the load parameter p, is used to calculate a suit- 
able approximation (predictor) for {Z} at a different 


value of p. This approximation is then chosen as an ini- 
tial estimate for { Z } in a corrective iterative scheme 
such as the Newton-Raphson technique. In each Newton- 
Raphson iteration the contact conditions are checked and 
updated. 

Computational Procedures to Determine 
Contact-Load Intensities, Contact Areas, 
and Energy Dissipated During Slip 

This section describes the contact algorithm used to 
assess the state of contact in the tire footprint and to 
implement the modified Coulomb friction law. Descrip- 
tions of the algorithms used to approximate the area of 
contact and the energy dissipated during slip are also 
included. Finally, the computational procedures used to 
determine the displacement, stress-resultant, and contact- 
load intensity solutions at each iteration and load step are 
outlined. 

Nonlinearities due to large displacements (moderate 
rotations) and the contact condition are combined into a 
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single iteration loop. Reference 23 advocates a two-level 
(nested) iteration scheme. For this two-level scheme, the 
inner iteration loop accounts for the contact conditions 
associated with the contact-load intensities, and the outer 
iteration loop uses the Newton-Raphson iteration 
scheme. Numerical experiments demonstrate that for 
frictionless contact problems the two-level iterative 
scheme requires more iterations than the single-level 
scheme utilized in the present study. (See ref. 29.) 

Contact Algorithm With Friction 


to 2. The dynamic friction-force load intensities are com- 
puted from the following equations: 


slip - COS ( V|/) ^vi M-dynamic 

(12a) 

slip = s i n (V)^nM- dynamic 

(12b) 


Incremental friction-force load intensities required to 
bring the friction-force load intensities back into compli- 
ance with the Coulomb friction law constraint are com- 
puted from 


Figure 3 schematically shows the contact algorithm 
developed for this investigation. Three possible contact 
states are possible for each contact node in this algo- 
rithm: (1) open or no contact, denoted by the contact flag 
set to 0; (2) stick contact where the contact node adheres 
to the contact surface, denoted by the contact flag set 
to 1; and (3) slip contact where the contact node slides on 
the contact surface, denoted by the contact flag set to 2. 
The algorithm is built on a two-level logical if-statement 
scheme. The first level interrogates the status of the pre- 
vious contact flags for each contact node. The second 
level interrogates the status of the current normal gap g w 
of nodes not previously in contact to determine whether 
or not contact has been established for that particular 
node; for nodes which were previously in contact, the 
second level interrogates both the sign of the normal 
Lagrange multiplier X w , which represents the normal 
contact load intensity, and the Coloumb friction law con- 
straints to determine the current nodal contact status. 

For nodes not previously in contact, stick contact is 
assumed whenever the normal gap g w is positive or 0, 
i.e., inequality equation (la) is satisfied, and no contact 
is assumed whenever g w is negative. For nodes previ- 
ously in contact, stick contact is assumed if both the 
static condition of compressive (or 0) normal tractions 
(eq. (lb)) and the static friction constraint (eq. (2a)) are 
satisfied; slip contact is assumed if the compressive nor- 
mal tractions condition is satisfied but the friction con- 
straint is not satisfied at a node. If the compressive 
normal tractions condition is not satisfied at a node, then 
that node is removed from the list of active contact 
nodes. 

When an open nodal contact condition is encoun- 
tered the contact flag is set to 0; if the previous contact 
condition was either stick or slip then the three load 
intensities {X M , \ v , X w } for that node are set to 0. When 
a stick nodal contact condition is encountered the contact 
flag is set to 1 , the number of boundary condition nodes 
is incremented, and incremental tangential nodal dis- 
placements Aw and Av are set to 0. When a slip nodal 
contact condition is encountered the contact flag is set 


AA, W — slip 03a) 

A* v = K iS]]p -X v (13b) 

The delta gaps associated with the sliding friction condi- 
tion are computed as follows 




u 


( { rlbar } - AX. h ) 
{dibar} 


A Sv 


({rlbar} - AA,„) 
{dlbar} 


(14a) 

(14b) 


where equations (14) are derived from equation (C44) in 
appendix C and the delta gaps replace the displacement 
solution. The vectors {rlbar} and {dlbar} are defined in 
appendix C. The incremental friction-force load intensi- 
ties (eqs. (13)) are used to update the contact solution 


— C — C r -C 

K u — -h 

(15a) 

“ C - C r~ C 

h v — X v + AX V 

(15b) 

and the delta gaps (eqs. (14)) are used to update the dis- 

placement solution 


U c = U c + e relax A g c u 

(16a) 

c c c 

V = * +£ relax A Sv. 

(16b) 


where £j. e i ax is a relaxation parameter with a magnitude 
0 < Ereiax < 1 . The function of 8^^ is to enhance the 
convergence characteristics of the contact algorithm. 
Index c ranges from 1 to the number of contact nodes. 

Determination of Contact Area 

To obtain an accurate approximation of the contact 
area of the Space Shuttle nose-gear tire under static load- 
ing conditions, an algorithm was developed which 
employs information from the initial geometry of the tire 
footprint and updates this geometry to accommodate 
footprint deformations associated with contact. Figure 4 
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illustrates major features of this contact area algorithm. 
The 9-noded element is conveniently subdivided into 
four quadrants with a node on the comers of each quad- 
rant. Coordinates of the quadrant comers are defined as 

i i i . 

x = Xq + u (17a) 


i i i , t . 

y = .yo + v ( 17b ) 

where x ^ and are the initial coordinates of the node 
and u i and v l are the tangential nodal displacement 
components. Index i ranges from 1 to the number of 
nodes in the element. As figure 4 indicates, the area of a 
typical quadrant is 


4 ^ e / a b b a be c b c a a c 

^quad = °- 5 (* >’ -•* y + x y -xy + X y -x y 

a c c a c d d c da ad . 

+ x y - x y + x y - x y + x y -x y ) 


(18) 


where the superscripts a , b , c, and d denote the four cor- 
ners of the quadrant, starting from the lower left and 
advancing counterclockwise around the quadrant. Nodal 
areas associated with contact nodes within a 9-noded ele- 
ment are 


^node 

A 2 

node 

71 node 

^node 

A 5 

74 node 

'''node 

A 7 

n node 
'''node 
71 node 


0.25A 


1 

quad 


0.25 (M qll ad + '''quad) 


0-25/l quad 

0.25 (/4 ^ad + ^quad) 

O- 25 ' 4 ^ 

°-25Mquad + <uad) 

°- 25A Ud 

0-25 ( q Ua d '''quad) 


0.25(A^ ua(1 4- Aq Uacl + A^ uat | + A^ uac j) ^ 


(19) 


When there is only partial contact within an element, 
i.e., only some nodes are in contact, then a variable 
weighting function is used to modify the nodal contact 
areas as follows: 


A 


i 

ct 


( 


0.5 


1.0 + 



k - 1 


node 




J 


( 20 ) 


where A* ct is the modified contact area associated with a 
node in contact, and ct k is defined as 


k 

ct 


1 1 if node* is in contact 
^0 if node* is not in contact 


( 21 ) 


where node* is a node on the periphery of the quadrant or 
quadrants defining the nodal areas. For the corner nodes 
of the 9-noded finite element shown in figure 4, n CJ is 3; 
for the midside nodes n ct is 5; and for the interior node 
n ct is 8. The peripheral nodes associated with comer node 
1 are (2, 8, 9); the peripheral nodes for midside node 2 
are (1, 3, 4, 8, 9); and the peripheral nodes for the interior 
node 9 are (1, 2, 3, 4, 5, 6, 7, 8). Finally, the nodal areas 
from each element are assembled into a global contact 
node array so that the nodal areas reflect the contribu- 
tions from all shared elements. 


Special Treatment for Tread Grooves 

To facilitate adequate modeling of the 3- 
circumferential groove tread pattern of the Space Shuttle 
nose-gear tire, modifications were made to the vector of 
initial normal gaps {gw 0 } and to the variable weighting 
function used to define the contact area. For nodes in the 
region of a tread groove, a positive number was added to 
gw 0 to prevent that node from contacting the surface. For 
this investigation the following modification was made: 


5 »'o + 3 ' 1 


( 22 ) 


where the superscript ranges over the node numbers 
associated with tread groove locations. Nodal areas for 
contact nodes adjacent to the center groove were com- 
puted directly from equations (19) and the variable 
weighting function (eq. (20)) was not used for those 
nodes. Nodal areas for contact nodes adjacent to either of 
the remaining tread grooves were computed with the 
variable weighting function (eq. (20)) modified to reflect 
fewer peripheral contact nodes. For example, suppose the 
right edge of the finite element shown in figure 4 coin- 
cides with the left edge of a tread groove such that nodes 
(3, 4, 5) are in the groove and never contact the surface. 
For this situation n ct for midside nodes (2, 6) is 3, and 
their peripheral nodes are (1, 8, 9) and (7, 8, 9), respec- 
tively. For the interior node (9) n ct is 5, and the periph- 
eral nodes are (1, 2, 6, 7, 8). 


Evaluation of Contact Forces From Load 
Intensities 

Solution of the governing discrete equations of 
the entire structure generates the nodal displacements, 
the stress-resultant parameters, and the values of the 
Lagrange multipliers at the contact nodes. For each 
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individual element in contact, the intensity of the 
contact-load intensities at a node {T u , T v , r„}_are^ equal 
to the values of the Lagrange multipliers { X v , X w } at 
the same node. Reference 30 states the contact-load 
intensities may be expressed as a function of the nodal 
forces { P l u , P l vy P l n } as follows: 


p[ 


T l 

u 


u 


• = f N l N J dQ< 

’J'J 

V 

J 

1 V 

. pi n . 

fl (e) 

. K . 


where N* are the shape functions used in approximating 
the Lagrange multipliers and the generalized displace- 
ments, and is the domain of the contact element. The 
indicated numerical integration in equation (23) was 
implemented with a Newton-Cotes quadrature formula- 
tion. The range of both i and j in equation (23) is from 1 
to the number of displacement nodes in the element. Ref- 
erences 29, 31, and 32 use this method in the analysis of 
frictionless contact problems. This method of extracting 
contact forces from the load intensities is consistent with 
the finite-element formulation of the problem. 

A possible drawback of this method is the fact that 
the shape functions N* are computed only for the un- 
deformed tire and may not reflect contact area variations 
associated with tire deformations at the tire-pavement 
interface. An alternative method of obtaining contact 
forces from the load intensities that alleviates this draw- 
back is expressed in the following equation: 


where w slip and v s i ip represent the tangential slip in the 
tire meridional ana circumferential directions, respec- 
tively. Index c ranges from 1 to the number of contact 
nodes and the index k ranges from 1 to the number of 
iterations at each load step. For a given contact node, the 
total slip associated with a load step is 

iter 

«siip = X C (26a) 

£ = 1 

iter 

- I 4p' (26b) 

k = 1 

The magnitude of these total slip values ip and i>£ lip 
should increase monotonically over the iterations that 
involve sliding contact. Energy dissipated through slip 
was computed at each iteration for each contact node as 

< 27a > 

Ki‘ p ’ = K, (27b) 

where the energy dissipation expressions ( t/slip u , ^/slip v ) 
should always be negative since the friction forces 

(T u Acp T v A ct ) oppose the slip. Energy expressions 

(eqs. (26)) are summed over the number of iterations to 
obtain an estimate of the total energy dissipated through 
sliding at each contact node and load step as follows: 




r i 

pC u 


< A C' T / 

K 

- = < 

(■ A c,T / > 

K, 




(24) 


^slip u - 

— c 

^sliPp = 


iter -ic 

A c, X (“slip 
*= 1 J 

r iter -,c 


A ct X ( V slip 7 'v>^ ) 


L * = 1 


(28a) 

(28b) 


where c ranges from 1 to the number of contact nodes. 
Reference 30 discusses other approaches for determining 
the contact forces. 

Evaluation of Energy Dissipated During Slip 

For contact problems involving friction it is impor- 
tant to account for any slip that occurs in the contact area. 
In this investigation the slip at each contact node is 
approximated for each iteration with the following 
equations: 

<Np’ = A Su \' ta < 2 *» 


Computational Procedures for Solving Tire 
Contact Problem 

The computational procedure used in the present 
study is summarized as follows: 

Preprocessing and Initial Calculation Phases 

Step 1 . Model tire geometry, evaluate stiffness coeffi- 
cients (ref. 31), and generate input data includ- 
ing transformation matrices. 

Step 2. Select estimates for the penalty parameters and 
assume the contact status at the selected con- 
tact nodes. 
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Step 3. Generate linear element arrays. 

Solution Phase 

Step 4. Solve inflation-pressure case without contact 
using Newton-Raphson iteration scheme. 

Step 5. Generate initial normal gap between the 
inflated tire configuration and flat surface at 
designated contact nodes. 

• Begin displacement incrementation loop. 

• Begin combined contact and Newton- 
Raphson iteration loop. 

Step 6. Generate nonlinear element arrays, eliminate 
the stress resultants and the Lagrange multipli- 
ers from the elemental equations, and assemble 
the left- and right-hand sides of the equations. 

Step 7. Solve equation (8) for the incremental 

displacements. 

Step 8. Update the response vector for displacements, 
stress resultants, and the Lagrange multipliers: 

{Z (r+,) } = {Z (r) } + { AZ (r) } (29) 

Step 9. Check the contact status and modify the con- 
tact conditions at each node as needed: 

if g < 0 and X w < 0, then the constraint is 
active 

if g w > 0 or A. vv > 0, then the contact constraint 
is inactive 

if T r = Jtjl + T 2 v )< 7'„)X static then the con- 
tact condition is stick 

if T r = J(T 2 u + T 2 v ) > r„|A static then the con- 
tact condition is slip 

When the previous normal contact constraint is 
inactive, proceed as follows: 

a. If the current normal contact constraint is 
also inactive, then continue. 

b. If the current normal contact constraint 
is active, add the active contact contri- 
bution to the list of nodes with stick con- 
tact, increment boundary conditions, and 
continue. 

When the previous contact constraint condition 
is stick or slip, proceed as follows; 

a. If the current contact constraint is 
stick, increment boundary conditions and 
continue. 


b. If the current contact condition is 
slip, define dynamic friction forces and 
continue. 

c. If the current normal contact constraint is 
inactive, remove the node from the list of 
active constraints, zero out contact forces, 
and continue. 

If any current contact constraints are different 
from the previous constraints or if there are any 
active slip constraints, return to step 6. 

Step 10. Check the convergence of the Newton- 
Raphson iterations: 

e = [{AZ}/{AZ}/{ zf { Z >] 1/2 <Tolerance (30) 
n 

where n is the total number of degrees of free- 
dom in the model and the tolerance is pre- 
scribed. If convergence is achieved, then 
compute the contact forces at each contact 
node with equation (24) or with the following: 



and continue. Otherwise return to step 6. 

Step 11. If the prescribed displacement is greater than 
the specified maximum displacement, then 
stop. Otherwise, add additional displacement 
and return to step 6. 

The mixed-formulation finite elements used in this 
study have nine displacement nodes and four stress- 
resultant nodes and are designated as M9-4 elements in 
table 1 . 

Comments on Mixed Models, Perturbed 
Lagrangian Formulation, and Computational 
Procedure 

The following comments regarding the mixed mod- 
els, the perturbed Lagrangian formulation, and the com- 
putational procedure used herein are in order: 

1. Nonlinear terms in the finite-element equations of the 
mixed model (eq. (7)) have a simpler form than those 
of the corresponding displacement model (eq. (9)). 

2. Equation (7) includes both the Lagrange multiplier 
approach and the penalty method as special cases, as 
follows: 
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a. As the penalty parameters E n and e t approach infin- 
ity, equation (7) reduces to those of the Lagrange 
multiplier approach. 

b. When Lagrangian multiplier terms are eliminated 
in equation (7), the resulting equations are identical 
to the penalty method. 

3. The perturbed Lagrangian formulation alleviates two 
of the drawbacks associated with the Lagrangian mul- 
tiplier approach and the penalty method, namely, 

a. The regularization term in the functional results 
in replacing one of the zero diagonal blocks in the 
discrete equations of the Lagrangian multiplier 
approach by the diagonal matrix [R]/e in 
equation (7). 

b. The contact condition is satisfied exactly by trans- 
forming the constrained problem to an uncon- 
strained problem through the introduction of 
Lagrangian multipliers (the terms 


J ^*u8 u | d£l> and J h w g w d£l 

in equation (6)) rather than approximately as in the 
penalty method. However, the presence of the regu- 
larization terms (the terms 





in equation (6)) results in replacing the contact con- 
dition by the perturbed condition: 


£[tf]U} + [0]'{*}-{s o } = 0 (32) 


4. An important consideration in the perturbed 
Lagrangian formulation and in any penalty formula- 
tion is the proper selection of the penalty parameters 
E n and E r With the foregoing mixed models the pen- 
alty parameters can be chosen independently of the 
element size without adversely affecting the perfor- 
mance of the model. Accuracy of the contact solution 
increases with increasing values of the penalty param- 
eter e n . However, for very large values of £„, the equa- 
tions become ill-conditioned and thus round off errors 
increase. For small values of e, the tire footprint 
becomes compliant, i.e., there is little or no slipping, 
and the calculated friction forces are artificially low. 
For very high values of E t the contact algorithm with 
friction may become ill-conditioned. 

5. The elemental arrays [F], [5], {G(X)}, 

{M(H, X)}, and {F} are evaluated numerically 


using a Gauss-Legendre formula. The arrays [Q], 
[/?], and {g 0 } are evaluated using a Newton-Cotes 
formula. In both cases the number of quadrature 
points used is the same as the number of displacement 
nodes in the element. This results in under-integrating 
the arrays [Q] and [R] and avoids the oscillatory 
behavior of the contact-load intensity that has been 
observed when the arrays are fully integrated. Note 
that the use of Newton-Cotes formula allows the 
contact-load intensities to be evaluated at the displace- 
ment nodes. Appendix E gives details on the shape 
functions for the Gauss-Legendre formula and the 
Newton-Cotes formula. 

Numerical Results 

Description of Finite- Element Models 

To develop the finite-element models used in the 
analysis of the Space Shuttle nose-gear tire, the cubic 
spline approximation of the outer meridional surface of 
the tire half cross section was discretized into 75 poten- 
tial node points as indicated in figure 5. From this popu- 
lation of possible nodes, a smaller number of nodes was 
chosen to approximate the tire cross section. To model 
the tire inflation response, a single strip of 30 finite ele- 
ments was used to approximate the complete tire cross 
section. This model employed 61 nodes to characterize 
the tire meridian, and 480 stress-resultant parameters and 
293 nonzero generalized displacement parameters were 
used to synthesize the tire inflation response. 

Finite-element models employed to analyze the con- 
tact behavior and friction characteristics of the Space 
Shuttle nose-gear tire used 41 node points in one half of a 
meridional cross section (81 nodes for the entire cross 
section), and these nodes are denoted as the circular sym- 
bols in figure 5. Nodes associated with the circumferen- 
tial tread grooves are also highlighted in figure 5. In the 
meridional direction, the tread area of the tire was mod- 
eled with the highest density of nodes and the sidewall 
and bead areas were modeled with progressively fewer 
nodes. This meridional node pattern was used for each of 
the two-dimensional finite-element tire models employed 
in this investigation. The circumference of the tire was 
divided into 240 possible node points, and a smaller 
number of nodes was chosen from that population to con- 
struct the tire finite-element models. To refine the mesh 
in specific areas such as the contact zone, a higher den- 
sity of nodes was chosen from the population in the spe- 
cific region of interest. 

Figure 6 shows a map of elements and node loca- 
tions for one of the models used to analyze the contact 
problem. Figure 6 shows an array of elements with 
40 elements in the meridional direction and 1 8 elements 
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in the circumferential direction. Numbers in the left and 
right margins of figure 6 denote the beginning and end- 
ing element numbers in specific rows. Black dots super- 
imposed over the square grid of elements denote the 
individual nodes of the finite-element model. Several 
individual elements are shaded and shown in an 
expanded scale to illustrate the node numbering sequence 
that is used to minimize the bandwidth for the finite- 
element models. The complicating factor here is that the 
elements in the circumferential direction are joined along 
the top and bottom edge. The numbering scheme that 
is illustrated in the example shown in figure 6 provides 
a minimum bandwidth for this tire model. For this spe- 
cific example the bandwidth is 1635. The six rows of ele- 
ments in the middle of the array, containing elements 1 
through 240, comprise the possible contact region for 
this model. Figure 6 also shows the location of the cir- 
cumferential tread grooves of the Space Shuttle nose- 
gear tire. 

Three different models were used in the analysis of 
the Space Shuttle nose-gear tire in contact with a flat 
plate. These models, denoted as model 1, model 2, and 
model 3, are depicted in figure 7. Each model employed 
480 elements in the region outside the contact zone 
(0 < 0.27c, 0 > 0.271). Model 1 included 240 elements in 
the contact region of the tire (-0.271 < 0 < 0.27t) for a 
total of 720 elements. (See fig. 6.) For model 1 there 
were 14 076 nonzero generalized displacement parame- 
ters, 23 040 stress-resultant parameters, and 3159 

contact-load intensity parameters. Model 2 used a refined 
mesh within the contact region with 480 contact elements 
and a total of 960 elements overall. For model 2 there 
were 18 776 nonzero generalized displacement parame- 
ters, 30 720 stress-resultant parameters, and 6075 

contact-load intensity parameters. Model 3 employed a 
more refined mesh in the contact zone with 960 con- 
tact elements and a total of 1440 elements. Model 3 

employed 28 152 generalized displacement parameters, 
46 080 stress-resultant parameters, and 1 1 907 contact 
load-intensity parameters. A single iteration for model 1 
required about 1 2 min on a Cray 2 computer, and a single 
iteration for model 3 required about 12 min on a Cray 
Y-MP computer. 

Convergence Characteristics and Performance of 

Contact-Friction Algorithm 

Relaxation parameter and penalty parameter 
effects. Table 2 shows the effect of variations in the 
relaxation parameter E^x on the convergence character- 
istics of the contact-friction algorithm. To study this 
effect four values of were evaluated over seven 

load steps. For a relaxation parameter value of 1.0 no 
convergence was obtained within 40 iterations beyond 


the second load step. It should be noted that the first load 
step, which computed the inflation solution, involved no 
contact and the second load step involved contact with no 
sliding. Load steps 3-7 involved some tire sliding for the 
conditions summarized in table 2. A relaxation parameter 
value of 0.75 required a total of 104 iterations to obtain 
converged solutions at the seven load steps; for a value of 
Ereiax °f 0-5, only 36 iterations were required to cover the 
same load range; and a relaxation parameter value of 
0.25 required a total of 149 iterations. It is obvious from 
these results that the choice of the relaxation parameter 
can have a profound effect on the convergence character- 
istics of the contact-friction algorithm. 

Data shown in table 3 illustrate the oscillating fric- 
tional load intensities that were observed when the relax- 
ation parameter was set to 0.75. In table 3 both the lateral 
friction load intensity and the drag friction load intensity 
are shown to change sign 13 times during the course 
of the 28 iterations required to obtain a converged solu- 
tion at the third load step. No oscillatory behavior was 
observed for the friction -force load intensities when the 
relaxation parameter was set to 0.5 or 0.25. 

Table 4 summarizes the effect of tangential penalty 
parameter variations on convergence characteristics of 
the contact- friction algorithm. To evaluate this effect 
four values of e t were studied over seven load steps. For 
a tangential penalty parameter value of 1.0 a total of 
31 iterations was required to obtain converged solutions 
over the range of load steps tested. At a tangential pen- 
alty parameter value of 10 12 the number of iterations 
required to obtain converged solutions over the same 
load range increased to 45. Thus, approximately a 
50 percent increase in iterations was required for con- 
verged solutions when z t was increased by 1 2 orders of 
magnitude. 

The physical significance of the tangential penalty 
parameter is that it represents the distributed tangential 
stiffness of the tire-contact interface and therefore serves 
a double role of penalty parameter and tire stiffness 
parameter in the contact- friction algorithm. When £ ; is 
set to a value of 1 .0, the footprint of the contact-friction 
algorithm is compliant; little or no slip occurs in the tire 
contact zone and the predicted friction-force load intensi- 
ties are very small. When z t is set to 10 12 the footprint of 
the contact-friction algorithm is extremely stiff, the tire 
must undergo slip to satisfy the Coulomb friction law 
constraint, and the calculated friction-force load intensi- 
ties may be too high. To illustrate this point more clearly, 
figure 8 shows the normalized tire slip dissipation ener- 
gies in the lateral and drag directions for two values of z r 
Figure 8 indicates that the amount of energy dissipated 
by the tire due to slipping in the footprint is strongly 
influenced by the choice of £,. Data presented in figure 8 
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also indicate that the contact-friction algorithm predicts 
much more energy dissipated from lateral slipping by the 
tire than from slipping in the drag direction for static 
loading cases studied in this investigation. 

Figure 9 presents the influence of normal penalty 
parameter magnitude on the accuracy of total strain 
energy and total contact force. Total strain energy is the 
integral of total strain energy density over the full extent 
of the tire carcass, while total contact force is the integral 
of normal contact-load intensity over the entire contact 
zone. Strain-energy ratio, denoted by the solid line, and 
contact-force ratio, denoted by the dashed line, are plot- 
ted as a function of the base 10 logarithm of the penalty 
parameter in figure 9. Results in figure 9 indicate that 
total calculated strain energy and total contact force are 
insensitive to variations in the normal penalty parameter 
over the range of 10 6 to 10 15 . It should be noted that for 
these tire-contact problems, total contact force is much 
more sensitive to a poor choice in the penalty parameter 
than total strain energy, as denoted by the percent errors 
on the two ordinate axes. 

Friction coefficient and bad step size effects . 
Table 5 summarizes the effect of static and dynamic fric- 
tion coefficient variations on convergence characteristics 
of the contact-friction algorithm. Three different static 
friction coefficients were evaluated to study these 
effects. |i stat i C = 0.3 represents a wet concrete runway 
condition; |i static = 0.6 represents a dry concrete runway 
condition; and p static = 1 0 represents a maximum fric- 
tion coefficient associated with aircraft tires in general. 
The dynamic friction coefficient for each friction state 
was taken to be 85 percent of the static value. These fric- 
tion coefficient values are consistent with a substantial 
aircraft tire friction database that has been acquired over 
a number of years. Reference 33 gives further details on 
experimental friction measurements for aircraft tires and 
empirical relationships for predicting aircraft tire friction 
responses. Data in table 5 indicate that about 30 percent 
more iterations were required to obtain converged solu- 
tions over the range of load steps tested for the low 
friction surface and the high friction surface than for 
the friction coefficients representative of the dry run- 
way condition. These results indicate that the contact- 
friction algorithm is robust enough to handle the range of 
friction coefficients normally experienced in aircraft tire 
applications. 

Table 6 summarizes the effects of varying load step 
size on the convergence characteristics of the contact- 
friction algorithm. This study evaluates performance of 
the algorithm over a normal tire deflection range from 0 
to 0.6 in. In the first case the tire deflection range was 
covered in 10 load steps; in the second case this deflec- 


tion range was covered in 5 load steps; and in the third 
case the deflection range was covered in 4 load steps. 
Case 1 required 52 iterations and computed a tire load of 
3998 lb for a tire deflection of 0.6 in. Case 2 required 
28 iterations and case 3 required 26 iterations to obtain a 
converged solution at a tire deflection of 0.6 in. The cal- 
culated normal tire load for case 2 was 0.6 percent less 
than the predicted load of case 1 and the calculated nor- 
mal load for case 3 was 0.9 percent less than that for 
case 1 . It appears that the contact-friction algorithm can 
operate over a range of step sizes without serious degra- 
dation in performance. 

Conclusions 

A computational procedure is presented for the solu- 
tion of frictional contact problems for aircraft tires. 
The Space Shuttle nose-gear tire was modeled using a 
two-dimensional laminated anisotropic shell theory 
which includes the effects of variation in material and 
geometric parameters, transverse-shear deformation, and 
geometric nonlinearities. Contact conditions were incor- 
porated into the formula by using a perturbed Lagrangian 
approach with the fundamental unknowns consisting 
of stress resultants, generalized displacements, and 
Lagrange multipliers associated with contact and friction 
conditions. The contact-friction algorithm is based on a 
modified Coulomb friction law. A modified two-field, 
mixed-variational principle was used to obtain elemental 
arrays. This modification consists of augmenting the 
functional of that principle by two terms: the Lagrange 
multiplier vector associated with normal and tangential 
node contact load intensities and a regularization term 
that is quadratic in the Lagrange multiplier vector. 

Shape functions used in approximating generalized 
displacements and Lagrange multipliers were selected to 
be the same and differ from those used to approximate 
stress resultants. Stress resultants and Lagrange multipli- 
ers were allowed to be discontinuous at the interelement 
boundaries. Nonlinearities due to large displacements, 
moderate rotations, and contact conditions were com- 
bined into the same iteration loop and were handled by 
using the Newton-Raphson iterative scheme. 

Numerical results are presented for the Space Shuttle 
nose-gear tire subjected to inflation pressure loads and 
combined inflation pressure and contact loads against a 
rigid flat plate. 

Results from this investigation lead to the following 
observations and conclusions: (1 ) the choice of the relax- 
ation parameter is critical to the performance of the 
contact-friction algorithm. If the parameter is too large, 
oscillating friction load intensities occur; if the param- 
eter is too small, many iterations are required for 
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convergence. (2) The tangential penalty parameter is a 
measure of tangential stiffness of the tire-contact inter- 
face and has a strong influence on energy dissipated by 
the tire due to slip. Normal contact-load intensity distri- 
bution and strain energy are insensitive to variations in 
the normal penalty parameter. (3) The contact-friction 


algorithm is robust enough to handle the range of friction 
coefficients associated with aircraft tire applications. 

NASA Langley Research Center 
Hampton, VA 23681-0001 
January 18, 1996 
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Appendix A 


Fundamental Equations of Shell Theory Used in Present Study 

Appendix A summarizes the fundamental equations of the Sanders-Budiansky type shell of revolution used in this 
study. Effects of laminated, anisotropic material response and transverse-shear deformation are included in these 
relationships. 

Strain-Displacement Relationships 


3 W 1 ( u ¥ 1,2 

5 5 R x 2 \R X 5 ) 2 y 

(Al) 


(A 2 ) 

= -S" + ( 3 . - t) + (^ ■ 3 - w lr 2 - M 

(A 3 ) 

K * = 3 A 

(A 4 ) 

d s r 1 

K 0 = ♦, + - r d eh 

(A 5 ) 

2K sO = + + 

(A6) 

2e *3 = ~^ +d s w + t 

(A 7 ) 

2e 93 = -^r + ; 3 0 M ' + < t , e 

(A8) 


where and £9 are extensional strains in the meridional and circumferential directions, 2e^0 is the in-plane shear strain, 
K s and K0 are bending strains in the meridional and circumferential directions, 2 k 5 0 is the twisting strain, 2 e s3 and 2£0 3 

3 3 

are transverse-shear strains, 3^ = — , 3 0 = , and <|> is the rotation around the normal to the shell, which is given by 


♦ = [-J a e M + ( a * + v) v 

Nonlinear terms that account for moderate rotations are underlined with dashes in equations (Al) to (A 3 ). 


(A 9 ) 


Constitutive Relations 

The shell is assumed to be made of a laminated, anisotropic, linearly elastic material. Every point of the shell is 
assumed to possess a single plane of elastic symmetry parallel to the middle surface. The relationships between the 
stress resultants and the strain measures of the shell are given by 
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where c ir /^, and d \j (i, j = 1, 2, 6) are shell stiffness coefficients. Nonorthotropic (anisotropic) terms are circled and 
dots indicate zero terms. 
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Appendix B 


Formulas for Elemental Arrays [F], [S], [G(X)] 9 [M{H y X)] 9 [F], [0], [F],and {g Q } 

Appendix B gives explicit forms of the elemental arrays [F], [5], [G(X)], [ M(H , X)], and [P]. The arrays are 
partitioned into blocks corresponding to the contributions from individual nodes or stress-resultant approximation func- 
tions. Expressions for the typical blocks are given in table Bl. The order of the strains is e s , Eg, 2e^, Kg 2k s q, 2£ v3 , 
and 2803 . The order of the nodal displacement parameters is w, v, w, (jjy, and <|> 0 . 


In table Bl, A and N are shape functions associated with stress-resultant components; N l and N J are shape 
functions for generalized displacements; 5 is the number of stress nodes in the element; m is the number of displacement 
nodes in the element; and £2^ is the element domain. The range of the indices k and / is 1 to s; the range of the indices i 


3 3 

and j is 1 to m. Dots in the matrices refer to zero terms; 3 = ^ ; and 3 fi = — . 

d s 


Quantities Ej and e 2 are defined in terms 


of 


U -N 


u 1 -> 

— d a w 

R 2 r 0 


(Bl) 

(B2) 
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Table Bl. Concluded 
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Appendix B also gives explicit forms of elemental arrays [Q ] , [R ] , and {g 0 }. Table B2 gives expressions of typi- 
cal partitions. In table B2, N l and N J are shape functions for Lagrange multipliers and generalized displacements and c 
is the number of nodal points in contact within the element. The range of the indices i and j is from 1 to c, and the range 
of the index i is from 1 to m\ <g w > is the unit ramp (or singularity) function defined as follows: 


n 

<g > 
ft vv 


si (8 W > 0 ) 

0 (*„<<>) 


(B3) 


where g w = -g w and n = 0 or 1. Vector {g 0 } contains components g u , g v , and <g w > for each contact node. Com- 
ponents g u and g are defined in figure 2(b). 


Table B2. Explicit Form of Typical Partitions of Arrays [Q ] , [/?] , and {# 0 } 


Array 

Number of partitions 
(or blocks) 

Formula for typical partition 

[Q 1 

m X c 

j N r N J <g w > 0 d£l 

\K] 

c xc 

- J N'N J <g w > () da 

n ( ° 

{s 0 } 

c 

J"' 

Sv 

<*»■> 

>d£l 
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Appendix C 

Derivation of Newton-Raphson Tangential Operator Equations 


Taylor’s Series Preliminary Discussion 

The governing differential equations for the tire contact problem are 


-/, */ • 
t 

(*) 

h 

(e) 

\ m ni xx 

(e) 

• 

s t • 

* » 

L 

< 

X 

• +• 

m n ,xh 

• 

> — < 

pP • 
. *0 . 


(*) 


= 0 


(Cl) 


In equation (Cl)// is the element flexibility matrix, dimensioned 58 x and partitioned into sxs submatrices; S{ is the 
element linear strain-displacement matrix, dimensioned 58 X m5 and partitioned into sxm submatrices; m ni is the ele- 
ment stress-displacement displacement matrix, dimensioned s3m3 X m3 and partitioned into smxm submatrices; and q c 
and r c are matrices associated with contact, dimensioned mxc3 and c x c3 and partitioned into mxc and cxc sub- 
matrices, respectively. Subvectors /z, x, and X representing stress, displacement, and contact-load intensity are dimen- 
sioned s , m, and c3, respectively. The normalized external load vector P is dimensioned m5 and the vector of initial gaps 
go is dimensioned c3. The penalty parameter is denoted as e and the loading parameter is denoted as p. Superscript t 
denotes transpose and superscript ( e ) indicates that the equations are developed at the element level. Equation (Cl) can 
be expressed as three functions by carrying out the indicated multiplication: 

f(h,x,X) = -f,h+s I x + ^m nl xx = 0 (C 2 ) 

g(h, Jt, X) = s^h + q c X + m nl xh - pP = 0 (C3) 


j(h,x,X) = q c x + ^jY + g 0 = 0 

Consider a pointy “reasonably” close to a root of a function fx). Taylor’s series expansion about jc 0 is 

f(x ) = /(A: 0 ) + (A--Jc 0 )/'(jc 0 ) + ^j(jr-A: 0 ) 2 /"(jc 0 )+ ... 


(C4) 


(C5) 


If fix) is set to zero, then x is a root of the right-hand side of the Taylor series expansion about jc 0 (taking linear terms 
only): 


o = f(x 0 ) + (x-x 0 )f'(x 0 ) 

(C 6 ) 

(x-x 0 )f'(x 0 ) = -f(x 0 ) 

(C7) 

/<* o) 

X = X Q 

f( x o) 

(C 8 ) 

, f(x o) 

X-X Q = 0 = ; 

/ (*„> 

(C9) 
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x represents an improved estimate of the root to replace jc 0 


(« + 1 ) n ~(n + 1 ) 

JC = X +8 


(CIO) 


Taylor’s series expansion of equation (Cl) about /i 0 , Jt 0 , and considering the linear terms only, leads to the following: 


f(h,x,X) = f{h Q ,x 0 ,k 0 ) + (h-h 0 )& + (x-x 0 )^ + (X-\o)^ (Cl 1) 

g(h, x, X) = g(h 0 ,x 0 , X Q ) + (h - h 0 )^ h + (x - X 0 )^ + ( X - vff (C 1 2) 


j(h, x, X) 


j(h Q , x# X Q ) + (h- h 0 )|j + (x - , 0 )| + (X - X 0 )^t 


(Cl 3) 


Case I — No Contact 

Consider first the case for no contact where equation (Cl) is reduced to 




s 


l 




2 m nl XX 


m nt xh 



= 0 


and rearranging the nonlinear terms 


-f l h + s,x + -m nl xx = 0 
s^h + m ni xh - pP = 0 


-f, [si+ l 2 m nl x J 


( s l + m nl x ) 



Applying the Taylor series expansion yields the following equations: 

- f i h 0 + s l x 0 + \ m nl x 0 x 0~ ( h - h o)f l + (* " x oX s l + m nl x Q> = 0 

S t h o + m n ,x 0 h 0 - pP + (h — /i 0 )(sj + m n ,x 0 ) + (x - x 0 )m nl h Q = 0 
Let 


(h - /2q) = Ah 


and 


(*-*o) = ^ 


(Cl 4) 

(Cl 5) 

(Cl 6) 

(Cl 7) 


(Cl 8) 
(Cl 9) 

(C20) 

(C21) 
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Equations (Cl 8) and (Cl 9) become 


"// ( S l + m nl X o) 

( S 'l + m n/ /l 0 


Ah 

Ax 


fl h O- s l x Q~2 m nl X O x O 
[ -s I h Q -m nl x 0 h 0 + pP j 


rh 

rx 


(C22) 


Equation (C22) defines the tangent operator for the Newton -Raphson iterative solution procedure for problems that do 
not involve contact. Looking at the first part of equation (C22) where 


1 


(C23) 


-fi A h + (s, + m nl x 0 )Ax = f,h 0 -s,x 0 --m nt x 0 x 0 = rh 
then multiplying by the inverse of the flexibility matrix 

- AA + f, \s, + m nl x 0 )Ax = h Q - f^s,x 0 - m nl x 0 x 0 = f' rh (C24) 

and then substituting into the second part of equation (C22) yields 

(•*/ + m nl x o)fl X< < s l + m nl x O> Ax + m nl h 0 Ax = ~ s\h Q - m n[ x Q h 0 + pP - (s, + m nl x 0 )f] X rh = rx - (s\ + m^x^f'rh 

(C25) 

Equation (C25) is the governing differential equation for problems that do not involve contact with the stresses elimi- 
nated. Equation (C25) can be rewritten by defining the following terms that are employed in the tire modeling code used 
in this investigation: 


[k 22 ] = m nl h Q 


-1, 


[* 12 bar] = (s t + m nl x 0 ) 

[*12*]' = ( s l + m nl x 0> 

(rhbar) = f ^ rh 


and equation (C25) now becomes 


([^ 12 5 ] [fc^bar] + [k 22 ])Ax = (rx) - (rhbar) 

Case II — Contact 

Next consider the case of contact by reintroducing the contact terms into equation (Cl 7) 


(C26) 

(C27) 

(C28) 

(C29) 

(C30) 




(s, + m nl x ) 


t 

<lc 


h 


• 

X 

> — < 

PP 

X 


So . 


= 0 


(C31) 


Applying the Taylor series expansion yields 


fl h O + s l x O + 2 m nl X 0 X 0-( h - h Q )f l + ( X - X o' ){S l + m nl X o) = 0 


(Cl 8) 
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s \h o + m nl x 0 h Q + q c X Q -pP + (h- h 0 )(s, + m n) x c Q ) + (x - x Q )m n/ h 0 + (X- X Q )q c = 0 

t r c t r c 

9c*o + 7^-0 -£ 0 + (*-*o)?c + (*- ^0)7 = 0 


Let 


(C32) 

(C33) 


(X-X Q ) = AX 

Equations (Cl 8), (C32), and (C33) become 


(C34) 


-f, 

(s l + m nl x 0 ) 

• 


Ah 


fl h Q ~ S I X Q - f»/V0 


rh 

(s, + m n t x 0 ) 

• 

4 C 

< 

Ax 

' = ’ 

-s , i h 0 - m nI x 0 h Q -4 c X 0 + pP 

► = < 

rx r 





AX 


y 


r ) 1 

• 

/ 

4 c 

r s 

E 

V. J 


- 4 C X 0 ~ 7^0 + #0 




The contact equation (the last part of eq. (C35)) can be written as 

(7 = - 4 c Ax - Wo - 7 *-0 + So 

and solving for AX yields 


(C35) 


(C36) 


^ = "(l) ’^-( 7) 'wo + [ j) 'so-K (C37) 

Substituting equation (C37) into the second part of equation (C35) redefines the governing differential equation to 
account for contact terms 


/I //■ 1 J ^ /f v—1 ^ 

(s\ + m n i x {) )f l C s, + m nl x 0 )Ax + m nl h 0 Ax-q c [jj q[ Ax = - s,h 0 - m n ,x 0 h Q + pP - (s, + m n{ x 0 )f, rh ~q c [jj (Sn~W 0 ) 


= rx-(s l + m nl x { 


(S0-^0) 


(C38) 


where the underlined terms are those associated with contact. Equation (C38) is identical to equation (8) in the main 
body of the paper. The following terms are associated with the contact solution in the tire modeling code used in this 
investigation: 


[dlbar] = (7) Xq ’c 


sy v— 1 

(rlbar) = f (g 0 - q[x 0 ) - X Q 


(C39) 

(C40) 


(rib) = 


7J (s 0 -q c x 0 ) 


(C41) 


and equation (C38) can be rewritten as 

([k ]2 s] t [k l2 bar] + [k 2 ~>]-q c [d\bar])Ax = {rx) - [fc^sj^rhbar) - <? c (rlb) (C42) 


25 



where the underlined terms are associated with contact. The stresses are recovered from the displacement solution 
through the following equation: 

( h ) = (rhbar) - [k n bar](x) (C43) 

and the contact-load intensities are recovered from the displacement solution with 

(X) = (rlbar) - f dlbar ] ( jc) (C44) 
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Appendix D 


Transformation of Elemental Arrays From Shell Coordinates to Global Cartesian Coordinates 

Transformation of the displacement components from the shell coordinates ( s , 9, x 3 ) to the global Cartesian coordi- 
nates (x, y, z) is expressed by the following equation: 

{X} (e) = m{X} (< ° (Dl) 

where [T] is a block-diagonal transformation whose submatrix at each node is given by 


[T] in) 

(5x5) 


> > > w > X X 

e s e Q e s x e e 0 0 

1 0 

4-0 

0 1 


(D2) 


> > * (e) 

where e s and eQ are tangential unit vectors in the s- and 9-directions, respectively, 0 is the null vector, and {X} 

and { X }*''* are generalized displacements in shell coordinates and global Cartesian coordinates, respectively. Note that 

rotation components <|) 5 and 4> 0 arc not transformed since the outer surface of the tire was chosen as the reference surface; 

therefore, § s and <{>9 do not appear in the contact conditions. 


Elemental matrices [5] and 


rdAfl 

_dX_ 


and the external load vector { P ) are transformed from the shell coordinates to 


the global Cartesian coordinates as follows: 


[si-Msim 


DM 

ax 




>V~dM- 

_ax. 


[T] 


(D3) 

(D4) 


{P}^[r]‘{P} 


(D5) 


Nonlinear vectors {G(X)} and {M(H, X)} are evaluated with displacement vector {X} expressed in terms of 
{X} at the end of each iteration cycle. 
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Appendix E 


Details of Shape Functions for M9-4 Finite-Element Model 

Appendix E presents the expressions for the shape functions used in the M9-4 finite-element model in terms of the 
local quadrilateral (or natural) coordinates £, and T|. In figure El the open circular symbols denote the nine nodes for the 
generalized displacements and Lagrange multipliers, and the filled symbols denote the interior nodes for the stress 
resultants. Figure El also lists the coordinates of the four comer nodes. 



(-L-1) 

Figure El. Schematic of the M9-4 finite element. 


Shape functions for the biquadratic approximations are given by the following system of equations: 

= ^T1($-1)(T1-1) 
n 2 = Ir,(i -^ 2 )(ti - 1) 

N 3 = ^T|(S+1)(T1-1) 

N 4 = ^+1)(1-T1 2 ) 

H 5 = ^Tl(^+ 1)01+ 1) 

V 6 = ^(l-4 2 )(ll + l) 

N 1 = 1)01 + 1) 

JV 8 = ^-1)(1-T1 2 ) 

N 9 = (1-5 2 )(1-t| 2 ) (El) 
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When the stress nodes are located at the center of each element quadrant, the corresponding bilinear shape functions 
for the stress resultants are given by 

N V = |(1 -24-2T1 +4^r|) 

N 2 = t(l+2£-2Tl-4$Ti) 

N 2, = |(1 + 2 ^ + 211 + 4 ^) 

^V 4 = 1(1 -24 + 2r| -4 ^ti) (E2) 

Elemental arrays [F], [S], (GfX)}, {M(H, X)}, and {F} are evaluated numerically using a Gauss-Legrendre 
formula, and the arrays [Q ] , [R ] , and {g 0 } are evaluated using a Newton-Cotes formula. Reference 34 describes these 
numerical quadrature formulas. 
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Table 1. Characteristics of Mixed Finite-Element Models Used in Numerical Studies 


Designation 

Number of 
displacement nodes 

Maximum number of 
Lagrange multipliers 

Number of 
parameters per 
stress resultant 

Number of 
quadrature points* 

M9-4 

3x3 

3x3 

2x2 

3x3 


♦All elemental arrays are evaluated using Gauss-Legendre quadrature formulas except for [2] , [R] , and {g^} , which are evaluated us- 
ing Newton-Cotes formulas. 


Table 2. Effect of Relaxation Parameter on Convergence of Contact-Friction Algorithm 
[e„ = 1.0 E+12; e, = 1 .5 E+06; |a static = 0.6; H dynamjc = 0.51] 


Step 

^relax 

Normal deflection, in. 

Normal load, lb 

Iterations 

1 

1.0 

0 

0 

6 

2 


.025 

86.63 

4 

3 


.05 

221.43 

(a) 

4 


.1 

515.07 

(a) 

5 


.2 

1256.76 

(a) 

6 


.3 

2112.22 

(a) 

7 


.33 

2315.18 

(a) 

Total 




10 a 

1 

0.75 

0 

0 

6 

2 


.025 

86.63 

4 

3 


.05 

221.43 

28 b 

4 


.1 

515.07 

15 b 

5 


.2 

1256.76 

22 b 

6 


.3 

2112.22 

25 b 

7 


.33 

2315.18 

26 b 

Total 




126 b 

i 

0.5 

0 

0 

6 

2 


.025 

86.63 

4 

3 


.05 

221.43 

4 

4 


.1 

515.07 

3 

5 


.2 

1256.76 

7 

6 


.3 

2112.22 

8 

7 


.33 

2315.18 

4 

Total 




36 

i 

0.25 

0 

0 

6 

2 


.025 

86.63 

4 

3 


.05 

221.43 

32 

4 


.1 

515.07 

20 

5 


.2 

1256.76 

27 

6 


.3 

2112.22 

28 

7 


.33 

2315.18 

32 

Total 




149 


a Did not converge in 40 iterations. 
b Oscillating friction forces. 


32 







































Table 3. Iteration History for Contact Node 36 Illustrating Oscillating Friction-Load Intensities 


[Ereiax = 0.75; £ n = 1.0 E+12; E t = 1 .5 E+06; |i static = 0.6; Adynamic = 0.51; Normal deflection = 0.05 in.] 


Iteration 

Contact flag 

K’ P si 

K> P si 

X v , psi 

Lateral friction 
slip energy, 
in-lb 

Drag friction 
slip energy, 
in-lb 

1 

1 

0.0000 E+00 

0.0000 E+00 

0.0000 E+00 

0.0000 E+00 

0.0000 E+00 

2 

2 

.1238 E+01 

.6313 E+00 

.2232 E-01 

-.7834 E-04 

-.9799 E-07 

3 

1 

.1412 E+01 

.6313 E+00 

.2232 E-01 

.0000 E+00 

.0000 E+00 

4 

2 

.1413 E+01 

-.7206 E+00 

-.2548 E-01 

-.4470 E-04 

-.5591 E-07 

5 

1 

.1326 E+01 

-.7206 E+00 

-.2548 E-01 

.0000 E+00 

.0000 E+00 

6 

2 

.1326 E+01 

.6759 E+00 

.2390 E-01 

-.2096 E-04 

-.2622 E-07 

7 

1 

.1370 E+01 

.6759 E+00 

.2390 E-01 

.0000 E+00 

.0000 E+00 

8 

2 

.1370 E+01 

-.6982 E+00 

-.2469 E-01 

-.1082 E-04 

-.1354 E-07 

9 

1 

.1348 E+01 

-.6982 E+00 

-.2469 E-01 

.0000 E+00 

.0000 E+00 

10 

2 

.1348 E+01 

.6871 E+00 

.2430 E-01 

-.5327 E-05 

-.6663 E-08 

11 

1 

.1359 E+01 

.6871 E+00 

.2430 E-01 

.0000 E+00 

.0000 E+00 

12 

2 

.1359 E+01 

-.6927 E+00 

-.2449 E-01 

-.2685 E-05 

-.3359 E-08 

13 

1 

.1353 E+01 

-.6927 E+00 

-.2449 E-01 

.0000 E+00 

.0000 E+00 

14 

2 

.1353 E+01 

.6899 E+00 

.2439 E-01 

-.1337 E-05 

-.1672 E-08 

15 

1 

.1356 E+01 

.6899 E+00 

.2439 E-01 

.0000 E+00 

.0000 E+00 

16 

2 

.1356 E+01 

-.6913 E+00 

-.2444 E-01 

-.6700 E-06 

-.8380 E-09 

17 

1 

.1354 E+01 

-.6913 E+00 

-.2444 E-01 

.0000 E+00 

.0000 E+00 

18 

2 

.1354 E+01 

.6902 E+00 

.2441 E-01 

-.3345 E-06 

-.4183 E-09 

19 

1 

.1355 E+01 

.6902 E+00 

.2441 E-01 

.0000 E+00 

.0000 E+00 

20 

2 

.1355 E+01 

-.6907 E+00 

-.2442 E-01 

-.1673 E-06 

.0000 E+00 

21 

1 

.1354 E+01 

-.6907 E+00 

-.2442 E-01 

.0000 E+00 

.0000 E+00 

22 

2 

.1354 E+01 

.6905 E+00 

.2442 E-01 

-.8365 E-07 

.0000 E+00 

23 

1 

.1355 E+01 

.6905 E+00 

.2442 E-01 

.0000 E+OO 

.0000 E+00 

24 

2 

.1355 E+01 

-.6906 E+00 

-.2442 E-01 

-.4183 E-07 

.0000 E+00 

25 

1 

.1354 E+01 

-.6906 E+00 

-.2442 E-01 

.0000 E+00 

.0000 E+00 

26 

2 

.1354 E+01 

.6905 E+00 

.2442 E-01 

-.2091 E-07 

.0000 E+00 

27 

1 

.1354 E+01 

.6905 E+00 

.2442 E-01 

.0000 E+00 

.0000 E+00 

28 

1 

.1354 E+01 

-.4183 E+00 

-.1479 E-01 

.0000 E+00 

.0000 E+00 

Total 





-.1655 E-03 

-.2066 E-06 
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Table 4. Effect of Tangential Penalty Parameter on Convergence of Contact-Friction Algorithm 


t^relax “ 0.5, E n — 1.0 E+12, Static “ 0.6, Adynamic ~ 0*5 1 ] 


Step 


Normal deflection, in. 

Normal load, lb 

Iterations 

1 

1.0 

0 

0 

6 

2 


.025 

86.63 

4 

3 


.05 

221.43 

3 

4 


.1 

515.07 

3 

5 


.2 

1256.76 

7 

6 


.3 

2112.22 

5 

7 


.33 

2315.18 

3 

Total 




31 

1 

1.5 E+03 

0 

0 

6 

2 


.025 

86.63 

4 

3 


.05 

221.43 

4 

4 


.1 

515.07 

3 

5 


.2 

1256.76 

7 

6 


.3 

2112.22 

5 

7 


.33 

2315.18 

5 

Total 




34 

i 

1.5 E+06 

0 

0 

6 

2 


.025 

86.63 

4 

3 


.05 

221.43 

4 

4 


.1 

515.07 

3 

5 


.2 

1256.76 

7 

6 


.3 

2112.22 

8 

7 


.33 

2315.18 

4 

Total 




36 

i 

1.0 E+12 

0 

0 

6 

2 


.025 

86.63 

6 

3 


.05 

221.43 

6 

4 


.1 

515.07 

5 

5 


.2 

1256.76 

7 

6 


.3 

2112.22 

7 

7 


.33 

2315.18 

8 

Total 




45 
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Table 5. Effect of Static and Dynamic Friction Coefficients on Convergence of Contact-Friction Algorithm 

[£r e i ax = 0.5; e n = 1.0 E+12; e, = 1.5 E+03] 


Step 

Friction coefficient 

Normal deflection, in. 

Normal load, lb 

Iterations 

1 

Mstatic "" 0*3 

0 

0 

6 

2 

Mdynamic = 0*26 

.025 

86.63 

7 

3 


.05 

221.43 

7 

4 


.1 

515.07 

5 

5 


.2 

1256.76 

8 

6 


.3 

2112.22 

8 

7 


.33 

2315.18 

4 

Total 




45 

1 

Mstatic “ 

0 

0 

6 

2 

M-dynamic — ^*5 ^ 

.025 

86.63 

4 

3 


.05 

221.43 

4 

4 


.1 

515.07 

3 

5 


.2 

1256.76 

7 

6 


.3 

2112.22 

5 

7 


.33 

2315.18 

5 

Total 




34 

1 

Mstatic — ^ *0 

0 

0 

6 

2 

Mdynamic = 0*®^ 

.025 

86.63 

6 

3 


.05 

221.43 

7 

4 


.1 

515.07 

5 

5 


.2 

1256.76 

8 

6 


.3 

2112.22 

8 

7 


.33 

2315.18 

4 

Total 




44 


35 


















Table 6. Effect of Load Step Size on Convergence of Contact-Friction Algorithm 
f^relax “ 0.5, = 1«5 E+03, — 1.0 E+12, Astatic = 0.6, Adynamic = 0*51] 


Step 

Normal deflection, in. 

Normal load, lb 

Iterations 

i 

0 

0 

6 

2 

.025 

86.63 

4 

3 

.05 

221.43 

4 

4 

.1 

515.07 

3 

5 

.2 

1256.76 

7 

6 

.3 

2112.22 

5 

7 

.33 

2315.18 

5 

8 

.4 

2806.26 

4 

9 

.5 

3435.70 

8 

10 

.6 

3998.32 

6 

Total 



52 

i 

0 

0 

6 

2 

.025 

86.63 

4 

3 

.2 

1255.77 

8 

4 

.4 

2779.85 

5 

5 

.6 

3974.51 

5 

Total 



28 

i 

0 

0 

6 

2 

.025 

86.63 

4 

3 

.3 

2110.36 

7 

4 

.6 

3962.25 

9 

Total 



26 
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Inflated profile 


Deformed under contact 


Contact region and 
normal contact 
load intensity 




(a) Normal gaps. 



(b) Transverse gaps. 


Figure 2. Schematic representation of gap terms associated with tire contact problems. 
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Figure 3. Schematic diagram of contact-friction algorithm logic. 
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Finite element: 




• Typical quadrant area: 


t = 0.5*(x 1 y 2 - x 2 y 1 + x 2 y 9 - x 9 y 2 4- x 9 y 1 - x*y 9 + x*y 9 - x 9 y* + x 9 y 8 - x 8 y 9 + x 8 y ! - x ! y 8 ) 


• Nodal areas: 

^ = 0 . 25 *^ 

A node = 0-25*(Aq Uad + A~ uad ) 

A node = 0-25*Aq Uad 

A node = °- 25 *( A 5uad+Aj uad ) 

A node = «-25*<ad 

A Le = °- 25 *( A quad+ A J U ad) 

A Le = °- 25 * A quad 

A node = °- 25 *( A ^ad+ A iuad) 

A node = 0-25*(Aq uad A quad A quad A quad) 


Figure 4. Major characteristics of contact surface area algorithm. 
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M9-4 finite element 


Sector 

Number of elements (M9-4) for — I 


Model 1 

Model 2 

Model 3 

-0.2n < 6 < 0.2k 

240 (40 x 6) 

480 (40 x 12) 

960 (40 x 24) 

0 < -0.2rc, 0 > 0.271 

480 (40 x 12) 

480 (40 x 12) 

480 (40 x 12) 

Total 

720 (40 x 18) 

960 (30 x 24) 

1440 (40 x 36) 


Figure 7. Finite-element models of Space Shuttle nose-gear tire used in present study. 
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Normalized 
tire slip 

dissipation energy 



Figure 8. Effect of tangential penalty parameter variation on slip energy dissipated by tire, = 0.5; 
£„ = 1 .0 E+12; |i s t a ti c = 0.6; Adynamic = 0.5 1 ; normal deflection = 0.3 in. 



Figure 9. Effect of magnitude of normal penalty parameter £„ on accuracy of total strain energy and contact 
force, pq = 300 psi. 
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